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Despite the importance of Si:P delta-doped wires for modern nanoelectronics, there are currently 
no computational models of electron transport in these devices. In this paper we present a nonequi¬ 
librium Green’s function model for electronic transport in a delta-doped wire, which is described by 
a tight-binding Hamiltonian matrix within a single-band effective-mass approximation. We use this 
transport model to calculate the current-voltage characteristics of a number of delta-doped wires, 
achieving good agreement with experiment. To motivate our transport model we have performed 
density-functional calculations for a variety of delta-doped wires, each with different donor configu¬ 
rations. These calculations also allow us to accurately define the electronic extent of a delta-doped 
wire, which we find to be at least 4.6 nm. 

PACS numbers: 73.22.-f, 73.63.-b, 73.63.Nm 


I. INTRODUCTION 

Phosphorus doping in silicon with atomic precision has 
led to a variety of new nanostructures on the silicon plat- 
forirP. Foremost among these are metallic wires that are 
“one atom tall and four atoms wide”^^. These wires are 
made using a 5-doping technique that combines scanning- 
probe lithography with molecular-beam epitaxjiSHll This 
technique achieves both high-density carrier concentra¬ 
tions and excellent two-dim ens ional (2D) conhnement of 
phosphorus atoms in silicorP*^. For example, 2D doping 
densities of one in four inside a (001) silicon monolayer 
{i.e. 0.25 ML) have previously been reportecP. These 
high donor densities result in spatially confined electron 
transport wh en an in-plane voltage bias is applied to the 
nanostructurdiSRil, Therefore, these systems have sim¬ 
ilar quantum mechanical properties to undoped silicon 
nanowires which are confined structurall y — " —I The elec¬ 
tronic properties of Si:P systems have applications for 
quantum computing and quantum communication tech- 
nologieJi^^^. In this paper, we examine some of these 
pro perties fo r a phosphorus in silicon (Si:P) 5-doped 
wir c ! ^ l ^°m ilMini using density-functional theory (DFT), 
effective-mass theory (EMT), the nonequilibrium Green’s 
function (NEGF) formalism, and tight-binding (TB) the¬ 
ory. 

The long spin-coherence times and large Bohr radius 
of the phosphorus donor electron in silicon make this 
material an interesting candidate for spin-based devices 
and other nanoscale electronicd^iltSil gj.p i5_(loped wires 
might be used a s the interconnects between stationary 
and flying qubitJ^^^^, l ow-re sistivity source-drain con¬ 
tacts for nanoelectronic^^MI^ qj- one-dimensional (ID) 
spin chains for confined magnon transporlPZHUI The 
modern 5-doped wire is a quasi-ID row of phosphorus 
atoms oriented in the [110] crystallographic direction, 
with a width of 1.54 nm in the lateral [lIO] direction 


which is equivalent to two dim er ro ws (2DR) on the re¬ 
constructed (001) silicon surfac^^lUl ju addition, because 
the placement of p hosphorus atoms on the (001) plane is 
indeterministicP^, a variety of donor configurations are 
possible within a realistic wire. 

Gomputational models of Si:P 5-doped wires are lim¬ 
ited to our density-functional model (Ref. [5T|) and the 
empirical TB model of Ref. [511 where the NEMO 3 D pack¬ 
age is used to describe the equilibrium electronic proper¬ 
ties of a Si:P 5-doped wir^. There are currently no com¬ 
putational models of electron transport in a 5-doped wire. 
However, electronic transport in an Si:P 5-doped layer 
has previously been investigated using a semi-classical 
modeP^. 

In this paper we calculate the electronic properties of 
a Si:P 5-doped wire for a variety of donor configurations. 
We investigate the effects of donor disorder and accu¬ 
rately define the electronic extent of a 5-doped wire. The 
current-voltage (I-V) characteristics of two 5-doped wires 
have recently been reported experimentallj®. Therefore, 
we expand on our density-functional model by using it 
to develop the first computational model for electronic 
transport in a 5-doped wire. We calculate the I-V char¬ 
acteristics for the two 5-doped wires which have been re¬ 
ported experimentally and investigate the change in these 
characteristics for wires with larger in-plane widths. This 
computational model has wide applicability because it 
scales easily to large system sizes, which are needed to 
simulate realistic devices. 


II. EQUILIBRIUM ELECTRONIC PROPERTIES 
OF <5-DOPED WIRES 

In this section we present the results of our density- 
functional calculations for the equilibrium electronic 
properties of a variety of 5-doped wires. The equilib- 
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rium properties are those for which the potential differ¬ 
ence between the source and drain contacts is equal to 
zero. Of particular interest are those properties relat¬ 
ing to the spatial confinement or electronic extent of the 
donor electrons perpendicular to the axis of the wires. 
We make an important distinction between two closely 
related concepts: the electronic width and the electronic 
extent of the wires. The electronic width is defined as a 
measure of spread for the probability density of the donor 
electrons perpendicular to the wire axis (for example, 
the full-width half maximum of this probability distribu¬ 
tion). This width is sometimes referred to as the “effec¬ 
tive electronic diameter” of the wirefl The electronic 
extent of the wires is defined as the minimum distance of 
lateral separation at which two wires do not affect each 
other’s equilibrium electronic properties. Therefore, the 
electronic extent is the minimum distance at which two 
wires must be separated so they behave exactly as they 
do in isolation. By contrast, the electronic width is the 
region within which it is most likely to find the donor 
electrons. 


A. Density-functional theory 

The SIESTA package was used to apply the Kohn- 
Sham self-consistent density-functional method in the 
generalized-gradient approximatiorP^^^^fil This method 
was used to calculate the equilibrium electronic proper¬ 
ties of the Si:P d-doped wire shown in Fig. This 
wire is made of a single row of phosphorus atoms that 
have been doped into one dimer row (IDR) of the re- 



[ 110 ] 

A 


[Il0]<—© 
[ 001 ] 


FIG. 1. (color online) (a) A 2D schematic of the single¬ 
row wire showing only part of the donor plane, and (b) a 
3D schematic of an orthorhombic supercell for the single-row 
wire, with silicon atoms (white spheres), phosphorus atoms 
(red spheres), the periodic boundaries of the supercell (blue 
lines), and a dotted line drawn inbetween where the two dimer 
rows would appear on a reconstructed (001) silicon surface. 
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FIG. 2. (color online) 2D schematics of double-row wires A-F 
(see labels) showing only part of the donor plane with silicon 
atoms (white spheres), phosphorus atoms (red spheres), the 
periodic boundaries of the supercell (blue lines), and a dotted 
line drawn inbetween where the two dimer rows would appear 
on a reconstructed (001) silicon surface. The periodic bound¬ 
aries are drawn such that there is always a donor atom at the 
origin of the supercells. 


constructed ( 001 ) silicon surface and will therefore be 
referred to as the single-row wire. The single-row wire 
represents the ID limit to scaling of the in-plane width 
of a phosphorus wire in silicorPS. A supercell for the 
single-row wire is shown in Fig. where the axis of the 
wire is oriented in the [110] direction. This supercell has 
dimensions of 0.77 nm x 5.40 nm x 5.46 nm, which cor¬ 
responds to at least 2.7 nm of bulk silicon “cladding” 
perpendicular to the wire axi^^. To avoid surface effects 
in the calculations, periodic boundary conditions are ap¬ 
plied to the supercells. In Fig.[l] the periodic boundaries 
are shown as blue lines. The length of the supercell in the 
[ 110 ] and [ 001 ] directions is therefore determined by the 
amount of bulk silicon cladding that is needed to isolate 
the phosphorus wire from its periodic images. 

We have also performed density-functional calculations 
on the Si:P d-doped wires shown in Fig. These wires 
have been doped into two dimer rows (2DRs) of the re¬ 
constructed ( 001 ) silicon surfac^^and will be referred to 
as double-row wires A-F (see labels in Fig. [^. Double¬ 
row wires A-F represent all donor configurations which 
can result from adsorption of PII 3 molecules onto a 2DR 
wide and 0.77 nm long region on the Si(OOl) surface, 
which is also periodic in the [110] directiorp3. For double¬ 
row wire B, we have previously found that by using 
2.7 nm of bulk silicon cladding perpendicular to the wire 
axis, the minima of the occupied bands in the resulting 
band structure are converged to within 5 meV^^. 

To reduce the periodicity of the donor atoms in the 
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[110] direction we would need to at least double the 
length of the supercell in this directiorP^. When the 
length of the supercell is doubled in one direction, the 
number of atoms in the supercell also doubles. This re¬ 
sults in an impractical eightfold increase in the compu¬ 
tation time of density-functional calculations. Therefore, 
the length of the supercells in the [110] direction is limited 
to 0.77 nm for all wires considered herein. For double¬ 
row wires A-F, the lengths of the supercell in the [110] 
and [001] directions are the same as that for the supercell 
shown in Fig. [^. However, the length of the supercell 
in the [110] direction is larger so that there is always at 
least 2.7 nm of bulk silicon cladding perpendicular to the 
wires. 



An optimized double-^ polarized basis set of localized 
atomic orbitals and norm-conserving Troullier-Martins 
pseudopotentials were used to variationally solve f or the 
ground-state electron density of these ^-doped wire^^MO] 
The lattice constant of bulk silicon was found to be 
5.4575 A using this basis set, which is in good agreement 
with the experimental value of 5.431 A^Sl. The exchange- 
correlation energy was calculated using the Perdew- 
Burke-Ernzerhof exchange-correlation functionaP^. To¬ 
tal energies were converged to within 10“"^ eV using a 
Methfessel-Paxton occupation function of order 5 with 
an electronic smearing of 0.026 eV. The mesh energy 
grid cut-off used was 300 Ry. A 6 x 1 x 1 Monkhorst- 
Pack /c-point grid was used to sample the Brillouin zone 
(BZ), which has previously been shown to give good re¬ 
sults for double-row wire #il. The atomic positions of 
the silicon atoms were not geometrically optimized after 
phosphorus substitution because their p osition s are not 
significantly affected by this substitutioiP^HH. Further 
details of our density-functional calculations are available 
in Refs. |3T1 [321 and |331 
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FIG. 3. (color online) The local density of states for the single¬ 
row wire (integrated between Fi and the equilibirum Fermi 
level and then line-averaged in the [001] direction) showing 
the probability density of the donor electrons (dark) with 
marked positions for phosphorus atoms (x), out-of-plane sili¬ 
con atoms along the wire axis (o), and in-plane silicon atoms 
(+)■ 


FIG. 4. The band structure for the single-row wire between 
r and |XoRT (solid lines) with the equilibrium Fermi level 
(dashed line) and the band structure for bulk silicon (shaded 
region). The GB minimum of bulk silicon has been set to en¬ 
ergy zero and the point Xort is at ^ in the [110] fc-space 
direction. For means of comparison, this band structure has 
been calculated using a supercell with the same dimensions 
as the supercells of double-row wires A, B, D, and E. 


B. Probability density of donor electrons 

The probability density for the donor electrons in the 
single-row wire is shown in Fig. [^ This probability den¬ 
sity has been calculated by integrating the local density of 
states (LDOS) between the conduction band (CB) edge 
and the equilibrium Fermi level of the single-row wire 
band structure shown in Fig. [^ Therefore, this prob¬ 
ability density is a mixture of the states with minima 
Fi, r 2 , and Ai (which are labeled in Fig. |^. We can as¬ 
sume the donor electrons occupy the bands in this energy 
range because the valence band of bulk silicon is fully 
occupied at equilibrium and so the only bands that are 
available for the donor electrons are those of the CB. Ex¬ 
perimental measurements of these systems are performed 
at T = 4.2 K and so we can also ignore the effects of elec- 
tro nic s mearing (due to thermal motion) in this calcula- 
tioiPni. 

The LDOS is calculated for the full supercell on 
a three-dimensional (3D) Cartesian grid, which is then 
line-averaged along the [001] direction (perpendicular to 
the donor plane) to make the probability density shown 
in Fig.[^ The probability distribution is normalized such 
that the integral of this probability density over the su¬ 
percell is equal to the number of donor electrons inside 
the supercell. In Fig. [^ a distance equivalent to two 
supercells is shown in the [110] direction. 

The donor electrons are partially delocalized along the 
axis of the phosphorus wire in Fig.[^ There is significant 
probability density not only on the phosphorus atoms 
but also the intervening silicon atoms both in- and out-of- 
plane. The majority of the probability density is localized 
to the atomic sites along the axis of the wire. Fig.[^shows 
there is strong spatial confinement of the donor electrons 
perpendicular to the axis of the wire. The probability 
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FIG. 5. (color online) The local density of states for double-row wires A-F (integrated between the Fi minimum and the 
equilibrium Fermi level and then line-averaged in the [001] direction) showing the probability density of the donor electrons 
(dark) with marked positions for phosphorus atoms (x), out-of-plane silicon atoms along the wire axes (o), and in-plane silicon 
atoms (-I-). 


density decays sharply with distance from the wire axis 
in the [110] direction. The majority of the probability 
density can seen to be localized to ±0.5 nm of the wire 
axis and, therefore, one may conclude that ~ 1 nm is 
a good approximation for the electronic extent of the 
wire. However, as we will demonstrate, this is not a valid 
approximation. Indeed, we will show that even ^ 2 nm 
is a poor approximation for the electronic extent of the 
single-row wire. 

Fig. [5] shows the probability densities for double-row 
wires A-F. These probability densities have been calcu¬ 
lated in the same way as Fig.|^ In each subfigure, there is 
significant overlap between the wavefunctions of the two 
rows of phosphorus atoms. The presence of this overlap 
is independent of whether the rows are staggered rela¬ 
tive to one another or aligned (compare the probability 
densities of double-row wires A-C with those of double¬ 
row wires D-F). If the overlap between the two rows of 
phosphorus atoms is large enough such that they behave 
as a single wire, then the electronic width of this wire is 
strongly dependent on the donor configuration. The elec¬ 
tronic width varies by as much as ~ 1 nm over double-row 
wires A-F and, therefore, we expect a similar variation 
in the width of a realistic wire. 

Overall, it is difficult to determine the electronic ex¬ 
tent of a (5-doped wire from these probability densities. 
Instead, we use the band structure of the wires and, in 
particular, energy splittings in the band minima to de¬ 
termine the electronic extent of a d-doped wire. 


C. Band structure of d-doped wires 

The band structure of bulk silicon calculated using a 
1280-atom orthorhombic (ORT) supercell is shown as the 
gray shaded region in Fig. There are two CB edges 
shown for bulk silicon, which are each doubly degenerate; 


one at F and the other at O.dGAloRT, where ATort = 
in the [110] A:-space direction. An explanation of 
the location of the bulk CB edges in this band structure 
is given in Appendix]^ The band structure of the single¬ 
row wire is shown as solid lines in Fig. The nuclear 
potential of the phosphorus atoms has moved the CB 
edges into the bulk bandgap region. The degeneracy of 
these CB edges or valleys has been broken, resulting in 
the appearance of four separate valleys. The minima of 
these four valleys are labeled as Fi, r 2 , Ai, and A 2 in 
Fig.0 

It is well-known for d-doped systems that an enhance¬ 
ment of the valley-orbit interaction due to spatial con¬ 
finement will break the six-fold degeneracy of the bulk 
silicon CB edge, resulting in a valley splitting (VS)P^®^. 
In Fig. the ID confinement caused by the donor po¬ 
tential has moved the CB edges into the bulk bandgap 
region and lifted their degeneracy by a VS. We report 
two different VSs for the single-row wire; one at F labeled 
Fi — r 2 and another at /c = 0.46 Vort labeled Ai — A 2 . 
The magnitude of the Fi — r 2 splitting is found to be 
109 meV, which is within 10 meV of the same VS re¬ 
ported for (5-doped layer^^. It has also been reported for 
(5-doped layers that CB valleys with higher curva ture are 
moved further into the bulk bandgap regiorP^^^ and this 
is shown for the single-row wire in Fig. In addition, 
we find the VS is larger for CB valleys with higher cur¬ 
vature. The magnitude of the Ai — A 2 splitting is found 
to be 21 meV, which is 19% of the Fi — r 2 splitting. 

The position of the Fermi level shows three of the four 
CB valleys to be occupied at equilibrium; the Fi, r 2 , and 
Ai bands. The Ai and A 2 minima at fc = 0.46Vort in 
the [110] fc-space direction are symmetrically equivalent 
to two CB edges ai k = — 0.46Vort- This is due to 
the symmetry of the path F —>■ Vqrt which is discussed 
in Appendix Therefore, for the single-row wire, we 
predict four conducting modes to be available for electron 
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FIG. 6 . The band structure for double-row wires A-F between F and IAort (solid lines) with the equilibrium Fermi level 
(dashed line) and the band structure for bulk silicon (shaded region). The CB minimum of bulk silicon has been set to energy 
zero and the point Xqrt is at in the [110] fc-space direction. For double-row wires C and F, the location of the bulk 

silicon CB minimum along the path F —> Aqrt has changed because the length of the supercell in the [lIO] direction is larger 
for these wires (see Appendix 0 


transport at low voltage biases. 

The band structures for double-row wires A-F are 
shown in Fig. The location of the bulk silicon CB 
edge for double-row wires C and F is at fc « O.SSATort 
not k « 0.46Xort, as discussed in Appendix [X} For the 
double-row wires, the larger nuclear potential of more 
phosphorus atoms moves the CB valleys further into the 
bulk bandgap region compared to the single-row wire. 
The minima of these valleys are again labeled as Fi, r 2 , 
Ai, and A 2 . For all double-row wires (except double¬ 
row wire D) there is a third CB valley minimum at the F 
point, labeled as Fa in Fig. The position of the Fermi 
level shows five CB valleys (including F 3 ) to be occu¬ 
pied at equilibrium for double-row wires A, C, E, and F. 
We expect there to be eight available conducting modes 
for the double-row wires, which is twice the number of 
available conducting modes as the single-row wire. How¬ 
ever, there are only seven available conducting modes for 
double-row wires A, C, E, and F and six for double-row 
wires B and D when the symmetry of the path F —>■ Aqrt 
is taken into account. 

The magnitudes of the Fi — F 2 and Ai — A 2 splittings 
are different for each of the double-row wires and, there¬ 
fore, these VSs must be dependent on donor configura¬ 
tion. It has previously been reported for d-doped layers 
that the magnitude of the Fi — F 2 splitting is dependent 
on the in-plane configuration of the donor atom^^. Fig.[^ 
shows the magnitude of the VS is affected by whether the 
two rows of phosphorus atoms are aligned (A, B, C) or 
staggered (D, E, F) relative to one another. The Fi — r 2 
and Ai — A 2 splittings are largest for double-row wires 


A and D. These wires represent the donor configurations 
for which the two rows of phosphorus atoms are later¬ 
ally separated by the smallest distance. The Fi — r 2 
and Ai — A 2 splittings decrease as the distance of lateral 
separation increases in Fig. 


D. Electronic extent of a (5-doped wire in the 
one-dimensional limit 


We expect there to be eight available conducting modes 
for double-row wires A-F. However, in the previous sub¬ 
section, we report seven (and six) available conducting 
modes for double-row wires A, C, E, and F (and double¬ 
row wires B and D). Therefore, the number of available 
conducting modes in these wires is reduced by something 
so far unaccounted for. If there were zero overlap between 
the wavefunctions of each row of phosphorus atoms, the 
band structure for the double-row wires would be identi¬ 
cal to the band structure of the single-row wire (except 
with each band being doubly degenerate). Therefore, we 
suggest F 3 is a degenerate pair of Fi and that this de¬ 
generacy has been lifted by an energy splitting which is 
proportional to the wavefunction overlap between the two 
rows of phosphorus atom^^. We label this energy split¬ 
ting as Fi-r 3 . In Fig. excluding double-row wires A 
and D, the ri-r 3 splitting decreases as the lateral sepa¬ 
ration of the two rows increases, i.e. as the wavefunction 
overlap between the two rows decreases so too does the 
ri-r 3 splitting. A similar energy splitting has prev iously 
been reported for two adjacent d-doped layerd^^M] 

In Fig. the ri-r 3 splitting is plotted versus the lat- 
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Lateral separation of single-row wires (nm) 


FIG. 7. (color online) The Fi—r 2 and Fi—Fs splittings versus 
the lateral separation of two single-row wires when they are 
aligned and staggered with respect to each other. The Fi — F 2 
splitting for an isolated single-row wire is also shown (dashed 
line). 


eral separation of the two rows of phosphorus atoms (be¬ 
tween 0.8 nm and 4.6 nm). We also plot the ri-r 2 split¬ 
ting versus the lateral separation of the two rows of phos¬ 
phorus atoms in this figure. These energy splittings have 
been calculated by using supercells that are larger than 
those used for double-row wires A-F^il. In Fig. the Fi- 
r 2 splitting tends towards the value of the VS calculated 
for the single-row wire as the lateral separation of the 
two rows is increased. By contrast, the ri-r 3 splitting 
tends towards zero as the lateral separation is increased. 
This suggests Fa is a degenerate pair of Fi and that the 
wavefunction overlap between the two rows of phospho¬ 
rus atoms is what breaks this degeneracy. 

At a lateral separation of 1.9 nm, two additional CB 
valley minima appear at fc « 0.46Aort (when the two 
rows of phosphorus atoms are aligned), which we inter¬ 
pret as degenerate pairs of Ai and A 2 that have also 
been lifted by an energy splitting. These splittings tend 
to zero as the lateral separation of the two rows is fur¬ 
ther increased. When the lateral separation is equal to 
3.5 nm, a fourth CB valley minimum appears at F which 
we suggest is the degenerate pair of F 2 . The energy split¬ 
ting of the F 2 degeneracy is equal to 60 meV at a lateral 
separation of 3.5 nm and decreases to 37 meV by 4.6 nm, 
i.e. it has not converged to zero at a lateral separation 
of 4.6 nm. 

We can now use this analysis of the Fi — Fs splitting 
(and other energy splittings in r 2 , Ai, and A 2 ) to ap¬ 
proximate the electronic extent of the single-row wire. 
When the lateral separation of the two rows of phospho¬ 
rus atoms is large, the two rows will each behave as single¬ 
row wires. The electronic extent of the single-row wire 
is then the lateral separation at which the energy split¬ 
tings in Fi, r 2 , Ai, and A 2 are equal to zero. The ri-r 3 
splitting and the energy splittings of the Ai and A 2 de¬ 


generacies are approximately equal to 6 meV at a lateral 
separation of 4.6 nm. This is equal to the uncertainty in 
Fi, r 2 , Ai and A 2 . Therefore, to within the uncertainty 
of our density-functional method, the energy splittings 
in Fi, Ai and A 2 are indistinguishable from zero and 
so the electronic extent of the single-row wire is approx¬ 
imately equal to 4.6 nm. However, because the energy 
splitting in the r 2 degeneracy has not converged to zero 
by 4.6 nm, this is but a lower bound on the electronic ex¬ 
tent of the single-row wire. Nonetheless, we expect this 
to be a good approximation because the occupancy of 
the Fi state is much greater than that of the r 2 state, 
as shown in Fig. For Vdoped wires that are separated 
by in-plane distances less than 4.6 nm, we predict a de¬ 
crease in the number of available conducting modes at 
low voltage biases. 

III. ELECTRONIC TRANSPORT PROPERTIES 
OF 5-DOPED WIRES 

In this section we use the results of our density- 
functional calculations to construct a computational 
model of electron transport in a Si:P (5-doped wire. We 
solve for the electronic transport properties of a (5-doped 
wire using the general NEGF approach of Datta and oth- 


A. The nonequilibrium Green’s function formalism 

The difficulty of solving for the eigenstates of a many- 
body system in the Schrodinger picture is avoided in the 
NEGF formalism by replacing the Hamiltonian operator 
by a Green’s function matrix. The transport properties of 
the system are then calculated from this Green’s function 
matrix. Our NEGF model describes a Si:P Vdoped wire 
using a TB Hamiltonian matrix, within a single-band 
effective-mass approximation, that is defined as 

N N 

H = ^ 

where e = —2Dt + U is the on-site energy, t = /2fha^ 
is the tunneling parameter, i and j are first-nearest- 
neighbor donor atoms, and N is the total number of 
donor atoms. U is an offset to the on-site energy due 
to a gate voltage applied to the wire (with G = 0 eV 
for a gate voltage equal to zero), D is the dimension 
of the device, fh is the effective mass of the donor elec¬ 
trons (discussed below), and a is the distance between 
two nearest-neighbor donor atoms. 

A Si:P (5-doped wire is divided into three parts: a 
source, a drain, and a channel region that separates the 
two. In general, Eq. [^describes the channel region only. 
However, in our NEGF model the source and drain con¬ 
tacts are described as semi-infinite extensions of the chan¬ 
nel and, therefore, Eq. [^is also used to describe the con- 
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tacts. The retarded Green’s function matrix is then de¬ 
fined as 


G(i;) = [(£; + fry)I-H-Ss-SD]-' (2) 

where E is energy, 77 is a positive infinitesimal real num¬ 
ber, and Ss and Sd are the self-energy matrices for the 
source and drain respectively. These are given by 

51s.d = '’■s,dSs.d'’'1d (3) 

where t is the coupling matrix between the contact and 
the channel, and g is the surface Green’s function for the 
contact. We calculate the surface Green’s functions using 
the iterative scheme of |Sancho et aLj which solves the 
accompanying Dyson equation to arbitrary precisioii^. 
The transmission function for the device is written as 


r(£;) = dxtr(rsGrDGt) (4) 


where d is the degeneracy of the single band (discussed 
below), and Fg and Fq are the broadening matrices for 
the source and the drain respectively. These are given by 


rs.D = i 


Ss.D - (Ss,d)^ 


(5) 


In the Landauer-Biittiker formalism, the current can then 
be calculated using the equation: 


I = I T (E) [fs {E) - /d (if)] dE ( 6 ) 

where q is the elementary charge of an electron, h is 
Planck’s constant and / (E) is the Fermi-Dirac distri¬ 
bution function for the contacts, defined as 


/s,D (E) 


1 

X -(- g(E — fJ,S,D)/kBT 


(7) 


In Eq. fee is Boltzmann’s constant, T is temperature 
(in Kelvin), and /is,D is the chemical potential of the 
source or drain which are given by 


MS = M 


VsD 

2 


( 8 ) 


and 


VsD 

Md = M-^ 


(9) 


with the equilibrium chemical potential of the wir^^ 
and Vg^D tbe source-drain bias voltage. For our calcula¬ 
tions, the source-drain bias voltage decays linearly across 
the channel region. 

The effective mass of the donor electrons is needed to 
fully specify the TB Hamiltonian in Eq. This effective 
mass can be calculated for the single-row wire from the 
occupied CB valleys in the band struture shown in Fig.|^ 
The band dispersion in the neighbo rhood of the CB valley 
minima is approximately parabolicP^^^. The curvature 


TABLE I. The values used for the free parameters in our 
NEGF model, mo is the free electron mass. 


rail mo 

mtimo 

a (A) 

Tm M (eV) 

d 

0.9163 

^ 0.1905 

7.718 

4.2 

0.12 

4 


Ref. EH 
Ref. [TT] 


/3 of this parabola is related to the effective mass m of 
the donor electrons through the equatioiP^ 


2fh 


I3k^ 


( 10 ) 


Therefore, the effective mass of the donor electrons can 
be calculated by fitting parabola to the CB valleys in 
Fig.S However, it has previously been shown for double¬ 
row wire B that the curvature of these bands do not 
change significantly from their bulk values and so we may 
use the transverse and longitudinal effective masses of 
bulk silicon without modificatioiP^. In Fig. the CB 
valleys at F have high curvature; they are described by 
the transverse effective mass rrit of bulk silicorl^. The 
CB valleys at jfcj « 0.46Xort have low curvature; they 
are described by the longitudinal effective mass mi of 
bulk silicorPi]. 

In Section HI] it was shown for the single-row wire that 
there were four conducting modes available for electronic 
transport at low voltage biases. Therefore, we describe 
the single-row wire by one conducting mode that is four¬ 
fold degenerate within a single-band effective-mass ap¬ 
proximation by setting d = 4 in Eq. |4| We set the effec¬ 
tive mass of this conducting mode equal to the average 
of the effective masses for the four occupied conducting 
modes. There are two occupied CB edges at F and two 
occupied CB edges at jfcj « 0.46Xort- The average of 
the effective masses is given by 


m = (2mt + 2^'^ /A ( 11 ) 


where there is a reduction of | in the longitudinal ef¬ 
fective mass (as discussed in Ref. [JT] and Appendix . 
The bulk silicon effective masses mi and mt, and the 
other parameters used in our NEGF model, are listed in 
Table |l| The equilibrium chemical potential /i is equal 
to the difference between Fi and the equilibrium Fermi 
level for the single-row wire in Fig.|4| The atomic spacing 
a is equal to the distance between two nearest-neighbor 
donor atoms in Fig. 2D wires with widths greater 
than that of the single-row wire are modeled as many 
adjacent single-row wires where the lateral separation of 
the donor atoms is equal to a and the tunneling parame¬ 
ter for nearest-neighbor atoms in adjacent wires is equal 
to t. Hardwall boundary conditions are applied in the 
dimension perpendicular to the axes of the wires. 



















B. I-V characteristics of 5-doped wires 


In this section we present I-V characteristics for two 
i5-doped wires that have recently been measured in ex- 
perimenlfSl. One of the wires has a width of 4.6 nm, 
which is equivalent to 6 dimer rows (6DRs) on the recon¬ 
structed (001) silicon surface. The other wire has a width 
of 1.5 nm, which is equivalent to 2DRs on the same sur¬ 
face. In addition, we present I-V characteristics for wires 
that have widths larger than 4.6 nm. For means of com¬ 
parison, all the wires have a channel region with a length 
of 47 nm (which is equal to the length spanned by the 
6 DR wire in experiment!ii|) . 

The I-V characteristics of the 6DR wire are shown in 
Fig.ii. The I-V curve for ?7 = 0 eV shows a linear re¬ 
sponse when a non-zero source-drain bias voltage is ap¬ 
plied to the 6DR wire. This linear response is character¬ 
istic of metallic conduction and is in good agreement with 
experimenlp]. The size of the current is approximately 
2 times greater than in experiment when [7 = 0 eV (see 
Fig. le of Ref. ITT]) . 

In Fig. 1^, the I-V curve for [/ = 0 eV shows a linear 
response when a bias voltage is applied to the 2DR wire. 
The size of the current is approximately 6 times greater 
than in experiment when [7 = 0 eV (see Fig. If of Ref. ED 
and is exactly half that of the 6DR wire. Therefore, 
when [7 = 0 eV there are double the number of available 
conducting modes in the 6DR wire than in the 2DR wire. 
The ratio of the two currents is also equal to two when 
U = —0.12 eV as shown in Fig. The doubling in 
the number of available conducting modes is confirmed 
by the transmission functions T{E) for the two wires in 
Fig. (compare the transmission function of the 6DR 
wire with that of the 2DR wire at E = fj, in Fig. [^). 

The gradient of these I-V curves is equal to the dif¬ 
ferential conductance G = dl/dVsu of the wires, which 
when divided by jh is equal to the number of conduct¬ 
ing modes that are available for electron transport. To 
model the application of a gate voltage to the wires, we 
add a non-zero offset energy [7 to the diagonal terms of 
the Hamiltonian matrijP^ describing the channel region 
(see Eq. . An applied gate voltage will change the oc¬ 
cupancy of the conducting modes in the wire and thereby 
change the conductance of the device. 

Beginning with the I-V curve for U = 0.13 eV in 
Fig. [^, the differential conductance of the 6DR wire 
increases as U is decreased to zero and then becomes 
increasingly negative. A negative offset energy in our 
NEGF model is equivalent to a positive gate voltage in 
experiment and, therefore, this trend is in good agree¬ 
ment with experimeniPl (where the differential conduc¬ 
tance increases as increasingly positive gate voltages are 
applied to the 6DR wire [see Fig. le of Ref. E]). A 
negative offset energy moves the unoccupied conduct¬ 
ing modes to lower energies, making them available to 
conduction electrons and thereby increasing the conduc¬ 
tance of the device. Alternatively, a positive offset en¬ 
ergy moves the occupied conducting modes to higher en- 



VsD (mV) 



VsD (mV) 

FIG. 8. (color online) Cnrrent versns source-drain bias volt¬ 
age for the (a) 6DR and (b) 2DR wire over a range of offset 
energies U. A negative offset energy is equivalent to a positive 
gate voltage in experiment. 


ergies, emptying these states of electrons and resulting 
in zero current because there are no longer any states 
available to conduction electrons (as shown in Fig. ||afor 
[7 = 0.13 eV). 

I-V curves for the 2DR wire are shown in Fig. for 
the same offset energies that are applied to the 6DR wire 
in Fig. [^. The relationship between these offset energies 
and the differential conductance of the 2DR wire is the 
same as that of the 6DR wire up to [7 = 0.11 eV. Beyond 
this point, the differential conductance for the 2DR wire 
increases when [7 = —0.30 eV and then remains constant 
as U is made more negative (as shown for [7 = —0.40 eV 
and [7 = —0.50 eV). Therefore, when [7 = —0.30 eV, the 
current in the 2DR wire has reached saturation because 
all of the conducting modes are occupied. By contrast, 
in the wider 6DR wire there remain conducting modes at 
higher energies that can be made available to conduction 
electrons by applying larger gate voltages. 

Finally, in both Figs. [8}r and 1^, there is a non-linear 
I-V response for U = 0.12 eV {i.e. for U = fi). When 
the offset energy is equal to the equilibrium chemical po¬ 
tential, the mean energy of the conduction electrons is 
equal to the energy of the CB edge. The application of 
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a source-drain voltage bias broadens the energy of the 
conduction electrons but it also increases the energy of 
the CB edge. And, therefore, it is possible for the rate 
of change in the energy of the CB edge to be greater 
than the rate of change in the energy broadening of the 
conduction electrons. If this is the case then although 
the current will increase as the source-drain bias volt¬ 
age is increased, the change in the current will decrease 
(as shown in Figs. and for U = 0.12 eV). This 
non-linearity is not the same as the non-linear response 
reported experimentally for the 2DR wire (see Fig. If of 
Ref. [TT|) . 

In general, there is good agreement between experi¬ 
ment and our results for the I-V characteristics of the 
6 DR wire. It is obvious that the results for the 6 DR wire 
are in better agreement with experiment than the 2DR 
wire. This is the case for both the size of the current 
and the change in the current versus source-drain bias 
voltage. Therefore, we conclude that the 6 DR wire is 
well-described by a ballistic model of electron transport, 
whereas the narrower 2DR wire is not. It is likely that 
a single-band effective-mass approximation is able to re¬ 
produce the low-temperature transport properties of the 
6 DR wire because the electron transport occurs at low 
voltage biases and, therefore, the higher energy modes 
of the silicon band structure are insignificant. This is 
promising for device simulation because our NEGF model 
scales easily to large system sizes. 

There are a number of approximations in our NEGF 
model that could explain the discrepancies between the¬ 
ory and experiment. For the 6 DR wire these include 
the value of the free parameters in Table |lj the bound¬ 
ary conditions, the average of the effective masses and 
the degeneracy of the single band in our effective-mass 
approximation. However, changing these properties will 
never reproduce the non-linear response reported exper¬ 
imentally for the 2DR wir^i^^. To reproduce this non¬ 
linearity, we suggest it is necessary to extend our NEGF 
model to include donor disorder and non-ballistic trans¬ 
port. In this paper we have chosen to present the simplest 
model of electron transport in a Si:P (5-doped wire and, 
therefore, leave the investigation of these approximations 
to the subject of future work. 

The transmission functions and equilibrium chemical 
potentials for wires with widths of 2DR, 4DR, 6 DR, SDR, 
and lODR are shown in Fig. Ep- In this figure, U is con¬ 
stant across all wires and is equal to -0.12 eV. The trans¬ 
mission functions and values of the equilibrium chemical 
potentials for each of these wires show that the number 
of available conducting modes increase as the width of 
the wires increase. Therefore, the I-V characteristics for 
the wider wires in Fig. have larger differential con¬ 
ductances than the narrower wires. The energy of the 
lowest conducting modes decrease as the width of the 
wires increase, as shown in Fig. Do. The spatial confine¬ 
ment of the donor electrons, perpendicular to the wire 
axis, decreases as the width of the i5-doped wire increases; 
thereby moving the conducting modes to lower energies 
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FIG. 9. (color online) (a) Current versus source-drain volt¬ 
age bias for wires with a variety of in-plane widths when 
U = —0.12 eV. (b) The transmission function T{E) versus 
energy E for each of these wires showing the position of the 
equilibrium chemical potential p when U = —0.12 eV. 


in a similar fashion to the eigenenergies of an infinite po¬ 
tential well. This decrease in spatial confinement also 
decreases the energy splittings between the conducting 
modes, which is shown in Fig. Ep as a narrowing of the 
steps in the transmission functions. Finally, it should 
be noted that the width of the steps in the transmission 
functions are related to the Fi — Fa splittin^^. Other¬ 
wise the transmission functions for the 2DR, 4DR, 6 DR, 
SDR, and lODR wires would vary only by a multiplicative 
constant. 


IV. CONCLUSIONS 

The band structure of a ^-doped wire has been cal¬ 
culated for a variety of donor configurations. We find a 
valley splitting at F, in agreement with previous density- 
functional calculations of d-doped layers. In addition, 
for (5-doped wires comprised of more than a single row 
of phosphorus atoms, we find another energy splitting 
at F. This energy splitting (Fi — F 3 ) is not caused by 
the valley-orbit interaction but by wavefunction overlap 
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between the adjacent rows of phosphorus atoms. The 
Ti — Fa splitting is then used to calculate the electronic 
extent of the single-row wire, which is found to be at 
least 4.6 nm. When two single-row wires are separated 
by in-plane distances less than 4.6 nm, the resulting en¬ 
ergy splittings reduce the number of available conducting 
modes and, therefore, conductance at low bias voltages. 

Furthermore, we present a nonequilibrium Green’s 
function model for electron transport in a Si:P (5-doped 
wire. We calculate the I-V characteristics of a variety 
of <5-doped wires which have different in-plane widths, 
achieving good agreement with experiment. These wires 
show a linear response to an applied bias voltage, which is 
characteristic of the metallic conduction that is observed 
in experiment. We also show that the conductance can 
be controlled by applying a gate voltage to the wires. In 
the future this NEGF model could be extended to include 
donor disorder and non-ballistic transport. 
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Appendix A: Band folding in (5-doped wires 

We analyze the band structure of the (5-doped wires 
with respect to the band structure of bulk silicon. It is 
well-known that silicon is an indirect bandgap semicon¬ 
ductor with a sixfold degenerate CB edge^. These six 
CB edges are located at /cq ~ 0-85^ in the first BZ of 
the face-centered cubic (FCC) Bravais lattice, one along 
each of the six (100) directions (where a is the lattice 
constant of bulk silicon). The band structure of bulk sil¬ 
icon calculated using a 2-atom FCC unit cell is shown in 
Fig. where Apcc is a point of high symmetry in the 
first BZ of the FCC unit cell and the path F —Apcc 
lies along one of the (100) directions in reciprocal spac^^. 
The CB edges are located at « O.SSApcc in Fig.[To|as 
|Afcc| = The band dispersion in the neighborhood 
of the CB edges is approximately paraboli(PiIi£l and, 
therefore, in reciprocal space these CB edges can be rep¬ 
resented in 3D by spheroidal su rfac es of constant energy 
centered at fcj^as shown in Fig. 111. The spheroidal sur¬ 
faces are anisotropic because the curvature of the band 
dispersion in silicon is anisotropic. 

Fig. also shows the band structure of bulk silicon 
calculated using an 8-atom simple cubic (SC) unit cell. 
From a comparison of the two band structures, we see 
the location of the CB edges and, therefore, spheroids in 
reciprocal space is dependent on the real-space unit cell 
that is used for the calculation, which is a result of band 
folding as discussed in our earlier work (see Appendices 
1 and 2 of Ref. l43|l . In Fig. 10 the BZ is folded about 
k = - due to a doubling in the length of the supercell 


to 5.46 A (for the SC unit cell). When the length of 
the unit cell is increased in one dimension, the length 
of the BZ in the equivalent reciprocal dimension is de¬ 
creased. The CB edges are thereby folded along their 
corresponding reciprocal space dimension towards the F 
poinlf^at k = (0,0,0). The CB edges have been folded 


from ko « 0.85^ to ko « 0.15f^ in Fig. 
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and this is 


in the [100] direction from 2.73 A (for the FCC unit cell) 


also shown in Fig. [TTJd by the translation oftne spheroids 
along each of the cardinal k axes towards k = (0,0,0). 

The band structure of bulk silicon calculated using a 
1280-atom orthorhombic (CRT) supercell (for example, 
see Fig. [^) is shown as the gray shaded region in Fig. [I] 
There are two CB edges at energy zero; one at F and 
the other at k « 0.16^. These are each doubly degen¬ 
erate and, therefore, represent four of the six CB edges 
of bulk silicon. The other two CB edges are located at 
k « —0.16^ and are not shown as they are symmetri¬ 
cally equivalent to those at fc « 0.16^. The locations of 
the CB edges in reciprocal space are dependent on the 
supercell that is used for the calculation and the result¬ 
ing band folding of the SC (and ultimately FCC) band 
structure. Therefore, the location of the CB edges for 
the 1280-atom CRT supercell can be predicted from the 
band structure calculated using the 8-atom SC unit cell 
and simple geometric arguments. 

To see how the SC band structure can be used to cal¬ 
culate the location of the CB edges for the 1280-atom 
CRT supercell, consider the simulation cell for an Si:P 
(5-doped layer. We use a 16-atom tetragonal (TFT) unit 
cell to represent (5-doped layers because the phosphorus 
atoms are doped in-plane at densities of 0.25 ML and the 
supercell needs to include at least four silicon atoms in 
the donor plane (so one of these atoms can be substi¬ 
tuted by a phosphorus atom for a doping density of one 

infourpa 

The TFT unit cell is rotated by 45° about the [001] axis 
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(a) 


(b) 


kz kz 




(c) (d) 


FIG. 11. (color online) A 3D representation of the six-fold degenerate CB minimum of bulk silicon showing the CB valleys as 
spheroidal surfaces of constant energy centered at fco for (a) a FCC unit cell and (b) the same representation for a SC unit cell, 
where the spheroids have been translated along each of the cardinal k axes towards k = (0, 0, 0). (c) A 2D representation of the 
six CB minima for a TET supercell, where the spheroids in (b) have been folded to the k^ky plane and (d) a ID representation 
of the six CB minima for an CRT supercell, where the ellipses in (c) have been folded onto the line k^ = ky. The width of 
these ellipses perpendicular to the line = ky go to zero in the limit as the length of the CRT supercell in the [110] direction 
tends to infinity. These surfaces have been colored as a guide to the eye only, there is no other meaning intended by the colors. 


compared to the 8-atom SC unit cell. This rotation does 
not affect the location of the CB edges in reciprocal space, 
only their relative position in the first BZ of the TET 
unit cell. The length of the TET unit cell in the z {i.e. 
[001]) direction is greater than that of the SC unit cell 
and this folds the CB edges in the kz direction towards T. 
If the length of the TET unit cell in the [001] direction 
is increased, it becomes a TET supercell. If this TET 
supercell is large enough to separate the 5-doped layer 
from its periodic images in the [001] direction, then the 
CB edges in the kz direction are folded to the T poinlP^!. 
This folding is shown as a “flattening” of the spheroids 
to the kxky plane in Fig. 11In this figure, there are CB 
edges located at T and k « 0.15^ in the [100] direction 
and, indeed, these have previously been reported for a 
TET supercelP^. 


The length of the simulation cell must also be large in 
the [110] direction for an Si:P (5-doped wire so the ID 
confinement of the donor electrons can be modeled accu¬ 
rately. In addition, the length of the simulation cell in 
the [110] direction will be different to the length of the 
simulation cell in the [001] direction (for atomistic sim¬ 
ulations) because of the different crystallographic sym¬ 
metries in each of these directions. Therefore, we must 
use an CRT supercell rather than a TET supercell to 
simulate the (5-doped wires. 

The CRT supercell is elongated in the [110] direc¬ 


tion compared to the TET supercell, which shortens the 
equivalent reciprocal space dimension of the first BZ (par¬ 
allel to the [110] /c-space direction and perpendicular to 
the line k^ = ky). The four CB edges at fc « 0.15— in 


Fig. 11: are thereby folded along the [110] fc-space direc¬ 


tion towards the line k^ = ky, as is shown in Fig. 111. 
Two of the six CB edges are located at fc « ^0.15^^ « 


0.11^ along the line k^ = ky in Fig. [ll[ l and another 
two at k « —0.11^ along the same line. There is also a 
^ contraction of the valleys for these four CB edges as a 
result of this folding and, ultimately, the 45° rotation of 
the CRT supercell (relative to the SC unit cell)P^. This 
contraction of the valleys causes an effective doubling 
in the curvature of the band dispersion and, therefore, a 
reduction of | in the corresponding effective masJ^. In 
addition, the CB edges are only approximately folded 
onto the line k^ = ky. This approximation is only ex¬ 
act in the limit as the length of the ORT supercell in 
the [110] direction tends to infinity. Therefore, the CB 
edge along the path T —>• Xqrt in Figure is located 


at fc « 0.46A1ort 


0.16 


271- 


rather than k 


0.11 


271- 


(where ATort = 2 ^^ fc-space direction). 

For double-row wires C and F, where the length of the 
ORT supercell in the [110] direction is larger, the same 
CB edge is located at fc « O.SSXqrt ~ 0-13^ in the 
[110] A:-space directiorPSl. 
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